(import (scheme base))
(import (scheme write))
(import (scheme file))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; These are the equations for the TAI to UTC offset. The values are: UTC
;; calendar year at start of interval; UTC calendar month at start of
;; interval; UTC calendar day at start of interval; offset-ref; drift; and
;; mjd-ref. The offset at the UTC time specified by mjd in each interval is
;; given by:
;;
;;   offset = offset-ref + drift * (mjd - mjd-ref).
;;
;; The order of the list is important; the earlier equations should be
;; first.

(define offset-equations
  '((1961 01 01 1.422818 0.0012960 37300)
    (1961 08 01 1.372818 0.0012960 37300)
    (1962 01 01 1.845858 0.0011232 37665)
    (1963 11 01 1.945858 0.0011232 37665)
    (1964 01 01 3.240130 0.0012960 38761)
    (1964 04 01 3.340130 0.0012960 38761)
    (1964 09 01 3.440130 0.0012960 38761)
    (1965 01 01 3.540130 0.0012960 38761)
    (1965 07 01 3.740130 0.0012960 38761)
    (1965 03 01 3.640130 0.0012960 38761)
    (1965 09 01 3.840130 0.0012960 38761)
    (1966 01 01 4.313170 0.0025920 39126)
    (1968 02 01 4.213170 0.0025920 39126)))

;; The last equation is valid until 1972-01-01 UTC, at which point the
;; offset if 10 seconds.

(define year-late   1972)
(define month-late  01)
(define day-late    01)
(define offset-late 10)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define year-start-list  (map (lambda (x) (list-ref x 0)) offset-equations))
(define month-start-list (map (lambda (x) (list-ref x 1)) offset-equations))
(define day-start-list   (map (lambda (x) (list-ref x 2)) offset-equations))
(define offset-ref-list  (map (lambda (x) (list-ref x 3)) offset-equations))
(define drift-list       (map (lambda (x) (list-ref x 4)) offset-equations))
(define mjd-ref-list     (map (lambda (x) (list-ref x 5)) offset-equations))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define year-end-list
  (append (cdr year-start-list) (list year-late)))
(define month-end-list
  (append (cdr month-start-list) (list month-late)))
(define day-end-list
  (append (cdr day-start-list) (list day-late)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define days-to-start-of-month
  `#(,(+ 0)
     ,(+ 0 31)
     ,(+ 0 31 28)
     ,(+ 0 31 28 31)
     ,(+ 0 31 28 31 30)
     ,(+ 0 31 28 31 30 31)
     ,(+ 0 31 28 31 30 31 30)
     ,(+ 0 31 28 31 30 31 30 31)
     ,(+ 0 31 28 31 30 31 30 31 31)
     ,(+ 0 31 28 31 30 31 30 31 31 30)
     ,(+ 0 31 28 31 30 31 30 31 31 30 31)
     ,(+ 0 31 28 31 30 31 30 31 31 30 31 30)
     ,(+ 0 31 28 31 30 31 30 31 31 30 31 30 31)))

(define (days-since-0000-00-00 year month day)
  ;; We assume the date is given in the proleptic Gregorian calendar.
  (let* ((leap-days-prior-to-year  (+ (floor-quotient (- year 1) 4)
                                      (- (floor-quotient (- year 1) 100))
                                      (floor-quotient (- year 1) 400)))
         (leap-year? (and (zero? (remainder year 4))
                          (or (not (zero? (remainder year 100)))
                              (zero? (remainder year 400))))))
    (+ day
       (vector-ref days-to-start-of-month (- month 1))
       (if (and leap-year? (> month 2)) 1 0)
       (* 365 year)
       leap-days-prior-to-year)))

(define (calendar-date-to-mjd year month day)
  ;; We assume the date is given in the proleptic Gregorian calendar. MJD 0
  ;; is 1858-11-17 00:00:00 UT.
  (- (days-since-0000-00-00 year month day)
     (days-since-0000-00-00 1858 11 17)))

(define mjd-start-list
  (map calendar-date-to-mjd year-start-list month-start-list day-start-list))
(define mjd-end-list
  (map calendar-date-to-mjd year-end-list month-end-list day-end-list))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(define offset-start-list
  (map
   (lambda (mjd-start offset-ref drift mjd-ref)
     (+ offset-ref (* drift (- mjd-start mjd-ref))))
   mjd-start-list
   offset-ref-list
   drift-list
   mjd-ref-list))

(define offset-end-list
  (map
   (lambda (mjd-end offset-ref drift mjd-ref)
     (+ offset-ref (* drift (- mjd-end mjd-ref))))
   mjd-end-list
   offset-ref-list
   drift-list
   mjd-ref-list))

(define mjd-epoch (calendar-date-to-mjd 1970 01 01))
(define mjd-late (calendar-date-to-mjd year-late month-late day-late))

(define tai-second-start-list
  (map
   (lambda (mjd-start offset-start)
     (+ (* (- mjd-start mjd-epoch) 24 60 60) offset-start))
   mjd-start-list
   offset-start-list))

(define tai-second-end-list
  (map
   (lambda (mjd-end offset-end)
     (+ (* (- mjd-end mjd-epoch) 24 60 60) offset-end))
   mjd-end-list
   offset-end-list))

(define next-tai-second-start-list
  (append (cdr tai-second-start-list)
          (list (+ (* (- mjd-late mjd-epoch) 24 60 60) offset-late))))

(define previous-tai-second-end-list
  (cons #f (reverse (cdr (reverse tai-second-end-list)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(when #t
  (with-output-to-file "lib/clock/tai-to-utc-offset-early.scm"
      (lambda ()

        (display ";; DO NOT EDIT!\n")
        (display ";; Generated automatically by generate-tai-to-utc-offset-early.scm.\n")

        (display "\n");

        (display "(define (tai-to-utc-offset-at-utc-day-early utc-day)\n")
        (display "  (cond\n")
        (for-each
         (lambda (mjd-start offset-ref drift mjd-ref)
           (display "   (")
           (write `(>= utc-day ,(- mjd-start mjd-epoch)))
           (display "\n")
           (display "    ")
           (write `(+ ,offset-ref (* ,drift (- utc-day ,(- mjd-ref mjd-epoch)))))
           (display ")\n"))
         (reverse mjd-start-list)
         (reverse offset-ref-list)
         (reverse drift-list)
         (reverse mjd-ref-list))
        (display "   (else\n")
        (display "    (error \"The TAI to UTC offset is not defined prior to 1961-01-01 UTC\"))))\n")

        (display "\n");

        (display "(define (tai-to-utc-offset-at-tai-second-early tai-second strict?)\n")
        (display "  (cond\n")
        (for-each
         (lambda (tai-second-start tai-second-end offset-start offset-end next-tai-second-start previous-tai-second-end)
           (when next-tai-second-start
             (when (< tai-second-end next-tai-second-start)
               (display "   (")
               (write `(and strict? (>= tai-second ,tai-second-end)))
               (display "\n")
               (display "    (error \"The TAI to UTC offset is not defined for this argument\" tai-second))\n"))
             (when (> tai-second-end next-tai-second-start)
               (display "   (")
               (write `(and strict? (>= tai-second ,next-tai-second-start)))
               (display "\n")
               (display "    (error \"The TAI to UTC offset is not uniquely defined for this argument\" tai-second))\n")))
           (let ((valid-tai-second-start
                  (if previous-tai-second-end
                    (max previous-tai-second-end tai-second-start)
                    tai-second-start)))
             (display "   (")
             (write `(>= tai-second ,valid-tai-second-start))
             (display "\n")
             (display "    ")
             (write `(+ ,offset-start
                      (* (- tai-second ,tai-second-start)
                         ,(/ (- offset-end offset-start)
                             (- tai-second-end tai-second-start)))))
             (display ")\n")))
         (reverse tai-second-start-list)
         (reverse tai-second-end-list)
         (reverse offset-start-list)
         (reverse offset-end-list)
         (reverse next-tai-second-start-list)
         (reverse previous-tai-second-end-list))
        (display "   (else\n")
        (display "    (error \"The TAI to UTC offset is not defined prior to 1961-01-01 UTC\"))))\n"))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

